display('loading data....')
PeakO=dlmread('scATAC_seq_RAd4.txt','\t',0,1);
X=dlmread('scRNA_seq_RAd4.txt','\t',0,1);
fid = fopen('scATAC_seq_RAd4.txt');
data = textscan(fid, '%s %*[^\n]');
fclose(fid);
PeakName=data{1,1};
fid = fopen('scRNA_seq_RAd4.txt');
data = textscan(fid, '%s %*[^\n]');
fclose(fid);
Symbol=data{1,1};
fid = fopen('RAd4_coupling_matrix.txt');
data = textscan(fid, '%s %s %f');
fclose(fid);
[d f1]=ismember(data{1,1},Symbol);
[d f2]=ismember(data{1,2},PeakName);
D=sparse(f1,f2,data{1,3},length(Symbol),length(PeakName));
%%%%%%%%%%%%%%
K=3;
tic
display('selecting parameters....')
[w2 h2]=nnmf(X,K);
[w1 h1]=nnmf(PeakO,K);
[lambda1 lambda2]=defaultpar_CoupledNMF(PeakO,w1,h1,X,w2,h2,D,0.5,0.01)
display('runing coupledNMF....')
W10=rand(size(PeakO,1),K);
H10=rand(K,size(PeakO,2));
W20=rand(size(X,1),K);
H20=rand(K,size(X,2));
W10=W10*diag(1./sqrt(sum(W10.^2)));
W20=W20*diag(1./sqrt(sum(W20.^2)));
[W1,H1,W2,H2,Score,detr]=NMF_cluster_sep2(PeakO,X,D,K,300,lambda1,lambda2,W10,H10,W20,H20,500,1);
%[W1,H1,W2,H2,Score,detr]=NMF_cluster_sep2_lap(PeakO,X,D,K,300,lambda1,lambda2,W10,H10,W20,H20);
[d S1]=max(H1);
[d S2]=max(H2);
t=toc
%%tSNE
display('seperate tSNE....')
Y=dlmread('tSNE_ATACseq.txt','\t');
Y1=dlmread('tSNE_RNAseq.txt','\t');
mymap = [0 0 0;134 0 179;251 133 51];
figure
scatter(Y(:,1),Y(:,2)',20,S1,'fill')
colormap(mymap/255)
figure
scatter(Y1(:,1),Y1(:,2)',20,S2,'fill')
colormap(mymap/255)
%%%
display('joint embeding....')
H1=diag(1./sqrt(sum(H1.^2,2)))*H1;
H2=diag(1./sqrt(sum(H2.^2,2)))*H2;
m=size(H1,2);
VH_f=[H1*size(H1,2) H2*size(H2,2)];
dlmwrite('H1_H2.txt',VH_f,'\t')
T22=tsne(VH_f');
[d S1]=max(H1);
[d S2]=max(H2);
SS=[S1 S2];
figure
scatter(T22(1:m,1),T22(1:m,2),20,SS(1:m),'s','fill')
hold on
scatter(T22(1+m:end,1),T22(1+m:end,2),20,SS(1+m:end),'fill')
colormap(mymap/255)